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Final  Report 

Low  Complexity  Track  Initialization  and  Fusion  for  Multi-modal  Sensor 

Networks 

Principal  Investigator:  Dr.  Qiang  Le 
Hampton  University 

The  objective  of  the  project  is  to  undertake  research  and  development  in  the  area  of  wireless  sensor 
networks  to  provide  an  effective  and  affordable  solution  in  support  increasing  requirements  for 
reconnaissance  and  surveillance.  A  wireless  sensor  network,  such  as  unattended  ground  sensors  (UGS) 
network  represents  a  set  of  inexpensive  sensors  of  various  types  (modes)  that  provide  an  effective  and 
affordable  solution  for  battlefield  surveillance.  For  example,  acoustic  sensors  could  obtain  bearing  angles 
of  targets,  and  proximity  sensors  report  the  presence  or  absence  of  targets  within  their  sensing  ranges.  The 
ultimate  goal  of  this  research  is  to  explore  fundamental  performance  bounds  that  determine  how  well  a 
sensor  network  can  resolve  and  localize  multiple  targets  as  a  function  of  the  operating  parameters  such  as 
sensor  density,  the  threshold  settings,  and  the  spacing  between  the  targets.  The  major  contribution  of  this 
project  is  the  development  of  the  particle-based  probability  density  (PHD)  filter  for  binary  measurements 
using  proximity  sensors. 

This  three  year  works  is  summarized  in  six  publications:  Aerospace  2010  [1]  and  Fusion  20 10 [2].  Army 
Science  2010[3],  Aerospace  201 1  [4],  SP1E  201 1  [5],  and  TAES  2012[6], 

The  Aerospace  2010  works  presents  the  maximum  likelihood  localization  (ML)  algorithm  for  multi-target 
localization  using  proximity-based  sensor  networks.  Proximity  sensors  simply  report  a  single  binary  value 
indicating  whether  or  not  a  target  is  near.  The  ML  approach  requires  a  hill  climbing  algorithm  to  find  the 
peak,  and  its  ability  to  find  the  global  peak  is  determined  by  the  initial  estimates  for  the  target  locations. 
This  paper  investigates  three  methods  to  initialize  the  ML  algorithm:  1 )  centroid  of  k-means  clustering,  2) 
centroid  of  clique  clustering,  and  3)  peak  in  the  1-target  likelihood  surface.  To  provide  a  performance 
bound  for  the  initialization  methods,  the  paper  also  considers  the  ground  truth  target  positions  as  initial 
estimates.  Simulations  compare  the  ability  of  these  methods  to  resolve  and  localize  two  targets.  The 
simulations  demonstrate  that  the  clique  clustering  technique  outperforms  k-means  clustering  and  is  nearly 
as  effective  as  the  1  -target  likelihood  peak  methods  at  a  fraction  of  the  computational  cost. 


The  Fusion  2010  works  documents  our  initial  attempt  to  determine  the  operating  parameters  of  a 
proximity  sensor  network  to  achieve  a  specified  probability  of  resolution  at  a  specified  target  separation 
rs.  The  probability  of  resolution  is  derived  by  considering  three  necessary  conditions  to  resolve  two 
targets.  These  conditions  consist  of  1)  node  proximity  to  the  targets,  2)  sufficient  number  of  detections 
from  nodes  close  to  the  target,  and  3)  the  existence  of  the  “between"  node  that  cannot  detect  either  target. 
These  conditions  lead  to  a  design  strategy  that  determines  the  necessary  sensor  density  and  threshold 
settings  to  achieve  the  desired  Pr  for  a  given  rs.  Simulations  demonstrate  that  at  the  designed  sensor 
density  and  threshold  values,  the  actual  percentage  of  targets  resolved  achieves  the  desirable  level  of 
resolution  for  moderate  to  large  target  separations. 

Army  Science  2010work  documents  our  initial  investigations  on  multiple  target  tracking  filters  in 
proximity  sensor  networks.  To  this  end,  we  implement  a  PHD  filter  and  compare  it  to  a  clairvoyant  PF 
filter.  During  the  measurement  update,  the  filters  do  not  associate  targets  with  the  binary  measurements 
due  to  the  fact  that  the  response  of  a  binary  sensor  is  intertwined  with  the  location  of  all  targets.  The 
simulations  that  compare  the  PHD  to  clairvoyant  processing  demonstrate  the  feasibility  of  the  PHD  to 
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both  estimate  the  number  of  targets  as  well  as  estimate  their  target  states.  One  useful  feature  of  the  PHD 
filter  is  that  it  considers  various  combinations  of  PHD  particles  when  performing  the  measurement 
update.  The  advantage  of  this  feature  was  demonstrated  via  the  simulations. 

Aerospace  201 1  work  further  documents  our  investigation  of  multiple  target  tracking  filters  in  proximity 
sensor  networks.  To  this  end,  we  implemented  a  particle -based  PHD  filter  and  compared  it  to 
ClusterTrack,  a  clustering-based  particle  filter.  The  simulations  that  compare  the  PHD  to  the  ClusterTrack 
method  demonstrate  the  advantages  of  PHD  method.  For  the  2-D  tracking  case,  the  localization  accuracy 
of  the  PHD  is  far  better  than  that  of  the  ClusterTrack.  Unlike  ClusterTrack,  the  PHD  is  able  to  estimate 
the  number  of  targets.  The  advantages  of  the  PHD  may  be  due  to  the  fact  that  it  is  based  principled 
Bayesian  approach;  whereas  the  ClusterTrack  is  more  ad-hoc  in  its  design. 

This  SP1E  2011  work  investigates  how  the  sensor  density,  sensing  range,  and  target  separation  affect  the 
ability  of  the  PHD  filter  to  estimate  the  number  of  targets  in  the  scene  and  to  localize  these  targets  (as 
measured  by  four  different  metrics).  Two  possible  measurement  models  are  considered.  The  disc  model 
assumes  target  detection  within  a  sensing  radius,  and  the  probabilistic  model  assumes  1/r2  propagation 
decay  of  the  source  signal  so  that  the  probability  of  detection  decreases  with  range  r.  The  simulations 
demonstrate  the  simplistic  disc  model  is  inadequate  for  the  PHD  filter  to  estimate  the  number  of  targets, 
and  the  filter  for  the  disc  model  exhibits  difficulty  to  localize  widely  separated  targets  for  low  sensor 
densities.  On  the  other  hand,  the  more  realistic  probabilistic  model  leads  to  a  PHD  filter  that  can 
accurately  estimate  the  number  and  locations  of  targets  even  for  small  target  separations. 


The  TAES  2012  work  documents  our  investigation  of  multiple  target  tracking  filters  in  proximity  sensor 
networks.  To  this  end,  we  developed  a  novel  formulation  of  the  PHD  filter  for  proximity  sensors.  We 
compared  the  PHD  against  two  clairvoyant  filters  that  assume  the  number  of  targets  is  known  a  priori.  For 
measuring  the  localization  errors,  we  considered  a  measure  that  attempts  to  separate  out  the  effects  of 
cardinality  errors.  The  simulations  show  the  advantages  of  the  crossover  feature  in  the  PHD. 


Two  possible  measurement  models  were  considered:  the  disc  and  probabilistic,  so  that  two  model 
matching  cases  and  two  model  mismatching  cases  were  considered.  While  the  disc  model  has  wide 
appeal,  the  probabilistic  model  better  represents  the  phenomenology  that  dictates  how  sensors  respond  to 
target  signals.  The  simulations  demonstrate  that  the  PHD  filter  would  be  unable  to  accurately  estimate  the 
number  of  targets  if  the  disc  model  was  assumed  in  the  measurement  update.  This  is  the  case  because  the 
hypothesis  that  explains  the  pattern  of  detection  reports  is  not  unique.  On  the  other  hand,  the  PHD  is 
accurate  when  using  the  probabilistic  model  in  the  measurement  update  because  of  the  uniqueness  of  the 
best  hypothesis  to  explain  the  pattern  of  detection  reports. 


The  unique  hypothesis  emerges  because  it  is  assumed  that  the  radiating  power  of  each  target  is  a  fixed  and 
known  quantity.  In  reality,  one  would  expect  that  this  power  fluctuates  from  target  to  target.  If  these 
fluctuations  are  unbounded,  the  probabilistic  model  becomes  completely  ambiguous,  and  it  would  be 
impossible  for  the  PHD  filter  to  accurately  estimate  the  number  of  objects  in  the  scene.  However,  we 
expect  that  these  fluctuations  are  bounded.  Future  work  will  investigate  how  allowing  for  bounded  spread 
of  target  power  levels  will  affect  the  performance  of  the  PHD  filter.  Other  future  work  will  investigate 
better  dynamical  models  to  improve  the  tracking  performance  including  target  motion  and  birth/death 
models,  and  other  cluster  methods  for  closely  spaced  targets. 
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Finally,  the  new  PHD  filter  in  this  paper  uses  a  profoundly  different  measurement  model  than  the 
traditional  PHD  filter.  When  the  sensor  measurements  are  statistically  independent  conditioned  on  the 
target  hypothesis  as  in  this  paper,  the  multisensor  version  of  the  PHD  is  straightforward.  However,  for  the 
traditional  PHD  filter,  the  “exact”  multisensor  version  is  very  cumbersome  to  implement.  Usually, 
approximations  are  incoiporated  in  the  multisensor  implementation  of  the  traditional  PHD.  These 
approximations  are  necessary  due  to  the  constraints  that  a  measurement  corresponds  either  to  one  unique 
target  or  clutter.  Actually,  the  measurements  are  usually  derived  from  raw  signals  through  signal 
processing,  e.g.,  direction  of  arrival  estimation.  As  an  alternative,  the  new  PHD  method  could  easily  be 
modified  to  consider  the  likelihood  of  raw  signals  given  the  target  hypotheses.  The  new  PHD  method 
does  suffer  from  approximations  in  the  particle  implementation  of  the  multiobject  target  density.  On  the 
other  hand,  approximations  of  the  PHD  for  sensor  fusion  also  lead  to  performance  loss.  Future  work  can 
investigate  whether  the  new  PHD  method  operating  over  raw  signals  from  multiple  sensors  has  utility 
compared  to  multiple  sensor  implementation  of  the  traditional  PHD  method. 


Since  FY2010  and  FY201 1  annual  reports  already  include  Aerospace  2010  [1]  and  Fusion  20 1 0[2] .  Army 
Science  20 1 0[3],  Aerospace  2011  [4],  SP1E  201 1  [5]  as  attachments.  In  this  final  report,  only  TAES  [12] 
is  included  as  attachment. 
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Probability  Hypothesis  Density-Based  Multitarget 
Tracking  for  Proximity  Sensor  Networks 

Qiang  Le  and  Lance  M.  Kaplan 


Abstract — This  paper  investigates  the  feasibility  of  a  mesh 
network  of  proximity  sensors  to  track  targets.  In  such  a  network, 
the  sensors  report  binary  detection/nondetection  measurements 
for  the  targets  within  proximity.  This  work  proposes  a  new 
probability  hypothesis  density  (PHD)  filter  and  its  particle 
implementation  for  multiple  target  tracking  in  a  proximity  sensor 
network.  The  performance  and  robustness  of  the  new  method  is 
evaluated  over  simulated  matching  and  mismatching  cases  for 
the  sensor  models.  The  simulations  demonstrate  the  utility  of 
the  PHD  filter  to  both  track  the  number  of  targets  and  their 
locations. 

Index  Terms — proximity  sensors,  PHD  filter,  multiple  target 
tracking. 

I.  Introduction 

Proximity  (or  binary)  sensors  are  simple,  low  power  devices 
that  report  whether  or  not  they  detect  a  target.  Typically, 
these  sensors  are  used  to  cue  more  sophisticated  sensors  that 
require  more  power  in  order  to  classify  and  track  the  targets. 
This  paper  investigates  the  feasibility  of  a  mesh  network  of 
proximity  sensors  to  track  targets.  In  such  a  network,  the  no¬ 
detection  report  is  as  valuable  as  a  detection  report.  Previous 
work  has  revealed  the  potential  of  target  localization  and 
tracking  for  a  single  target  [1],  [2],  [3],  or  for  a  ideal  disc 
sensing  model  [4],  [5],  [6],  [7],  [8]  where  [7],  [8]  does  consider 
multiple  target  tracking.  There  are  two  models  to  describe 
the  behavior  of  the  proximity  sensor  that  drive  localization 
algorithms.  Many  approaches  are  based  on  the  disc  model 
that  assumes  a  sensor  detects  a  target  when  and  only  when  the 
target’s  distance  to  the  sensor  is  below  a  sensing  radius  [4],  [5], 
[6],  [7],  The  disc  model  ignores  physical  phenomenology  such 
as  the  additive  effects  of  the  energy  radiating  from  multiple 
targets  will  extend  the  sensing  range.  Other  approaches  use 
a  probabilistic  model  for  the  sensors  that  is  based  on  a 
physical  energy  loss  as  the  energy  is  radiated  from  each  of 
the  targets  [1],  [2], 

For  multiple-target  tracking,  Singh,  et  al. ,  proposed  Clus- 
terTrack ,  which  is  a  clustering-based  particle  filter  approach 
to  track  targets  based  on  the  disc  model  [7],  [8].  While 
ClusterTrack  can  track  targets  in  real-time,  the  method  can 
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only  localize  all  targets  when  operating  in  a  batch  mode 
where  measurements  from  later  points  in  time  are  employed 
to  localize  tracks  at  a  given  time.  Furthermore,  it  assumes 
the  number  of  targets  in  the  scene  does  not  change  over  the 
batch  processing  interval.  He,  et  al.  proposed  a  multitarget 
tracking  approach  using  a  probabilistic  model  [9],  However, 
the  method  only  tracks  the  presence  of  targets  within  cells  of  a 
grid.  On  the  other  hand,  the  focus  of  our  work  is  to  localize  and 
track  multiple  targets  over  continuous  states  including  target 
location.  Our  recent  work  presents  a  maximum  likelihood 
(ML)  localization  algorithm  for  multitarget  localization  when 
the  number  of  targets  is  known  [10].  This  work  extends  [10] 
by  considering  the  more  general  application  of  target  tracking 
when  the  number  of  targets  is  unknown.  To  this  end,  the  likeli¬ 
hood  function  is  integrated  into  a  particle  filter  implementation 
of  a  modified  version  of  the  probability  hypothesis  density 
(PHD)  filter  as  described  in  [11], 

Traditional  particle  filter  approaches  in  multiple  target  track¬ 
ing  must  employ  measurement-to-track  association  techniques 
to  ensure  as  much  as  possible  that  target  state  estimates  are 
not  contaminated  by  wrong  measurements  [12].  In  contrast, 
the  PHD  filter  introduced  in  [11]  is  able  to  “sweep  the 
association  issue  under  the  rug”  because  the  PHD  is  actually 
tracking  the  density  of  target  presence,  i.e.,  PHD,  and  not 
actual  point  targets.  The  estimated  number  of  targets  in  the 
scene  is  simply  the  integral  of  the  PHD  over  the  target 
states.  However,  extraction  of  estimates  of  the  target  locations 
requires  a  clustering  method.  The  issues  of  how  to  transform 
the  PHD  into  point  target  estimates  are  discussed  in  [13]. 

The  applications  of  the  PHD  filter  have  been  seen  in  sonar, 
visual,  and  radar  tracking  [13],  [14],  [15],  [16],  [17],  [18].  For 
example,  [14],  [15]  implemented  the  PHD  filter  to  identify 
the  underwater  obstacles  that  would  need  to  be  avoided  in 
navigation  using  forward-scan  sonar  images.  In  [16],  the  PHD 
filter  was  employed  to  track  a  random  number  of  pedestrians 
in  image  sequences  and  derive  their  location  sequences.  In 
[17],  the  PHD  filter  was  applied  to  target  tracking  using 
both  range  and  Doppler  measurements.  More  generic  PHD 
filters  are  described  in  [13],  [18]  that  incorporate  Cartesian 
target  position  measurements  directly  or  through  bearing  and 
range  measurements.  To  the  best  of  our  knowledge,  the  PHD 
filter  has  not  yet  been  used  in  proximity-based  multiple  target 
tracking. 

These  PHD  filters  assume  that  the  measurements  are  in¬ 
dependent  and  associated  to  one  target.  As  a  result,  the 
measurement  update  of  the  PHD  has  a  simple  analytical  form. 
This  paper  takes  a  more  general  view  that  the  PHD  filter  is  a 
sequential  Bayesian  tracker  that  is  only  following  the  evolution 
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of  the  first  order  multitarget  moment  of  the  multiobject  density 
no  matter  the  dynamical  and  measurement  models  that  are 
incorporated.  A  motivation  for  the  PHD  is  that  following  the 
multiobject  density  itself  is  too  cumbersome  for  practical  im¬ 
plementation,  e.g.,  via  particle  filtering.  In  this  paper,  we  apply 
PHD  filtering  for  a  mesh  network  of  proximity  sensors  where 
the  measurements  do  not  directly  associate  to  a  single  target. 
As  is  shown  in  the  next  section,  the  response  of  the  proximity 
sensor  depends  on  its  location  relative  to  the  location  of  all 
targets  within  the  scene.  As  a  result,  the  measurements,  i.e.,  the 
binary  outputs  of  the  proximity  sensors,  are  intertwined  with 
all  targets.  Therefore,  tracking  with  binary  sensors  requires  a 
modification  of  the  likelihood  (measurement)  update  equation 
in  the  original  PHD  formulation.  In  this  paper,  we  derive  the 
modified  measurement  update  equation  by  using  the  same 
Bayesian  update  method  as  in  [11],  with  the  exception  that 
the  form  for  the  multisensor/target  likelihood  is  not  separable. 
Unfortunately,  the  new  likelihood  update  method  does  not 
admit  a  closed  form  analytical  solution.  Nevertheless,  the 
numerical  computation  via  particle  filtering  is  straightforward. 

This  paper  is  an  extension  of  a  series  of  our  earlier  confer¬ 
ence  papers  that  initially  presented  our  ideas  for  PHD-based 
multitarget  tracking  via  a  network  of  proximity  sensors.  The 
initial  method  was  introduced  in  [19],  and  the  method  was 
compared  with  ClusterTrack  in  [20].  The  ClusterTrack  particle 
filter  is  to  the  best  of  our  knowledge  the  only  other  method  to 
track  multiple  targets  over  a  continuous  space  via  proximity 
sensors.  An  initial  version  of  the  ClusterTrack  particle  filter 
was  proposed  for  the  disc  sensor  model  and  tested  on  1  -D 
data  in  [7],  In  [20],  we  demonstrated  that  unlike  that  version  of 
ClusterTrack,  the  PHD  method  was  able  to  accurately  estimate 
the  number  of  targets,  and  the  PHD  methods  were  able  to 
better  track  the  targets.  In  [21],  we  investigated  the  perfor¬ 
mance  of  the  PHD  method  as  a  function  of  sensor  density 
and  target  separation  of  two  targets.  That  work  demonstrated 
that  the  PHD  method  using  the  probabilistic  sensing  model 
is  able  to  estimate  the  number  of  targets  and  localize  them 
no  matter  the  target  separation  for  sufficient  sensor  density. 
This  paper  expands  upon  [19],  [20],  [21]  by  presenting  a 
more  thorough  description  of  the  new  PHD  filter  and  its 
particle  implementation.  Furthermore,  this  paper  includes  a 
new  investigation  of  the  performance  of  the  PHD  method  for 
both  model  matching  and  mismatching  conditions.  We  also 
compare  the  PHD  methods  against  the  more  recent  version 
of  ClusterTrack  [8],  Furthermore,  this  paper  investigates  the 
tradeoff  between  the  relative  weighting  on  the  reliance  on  prior 
measurements  and  current  measurements  through  the  use  of 
“innovative”  particles. 

The  paper  is  as  organized  as  follows.  Section  II  provides 
details  about  the  proximity  sensor  model.  Then,  Section  III 
provides  a  basic  review  of  finite  point  processes  and  presents 
the  new  PHD  filter  for  the  proximity  sensor  network  and  its 
particle  filter  implementation.  This  section  also  discusses  a 
clairvoyant  A' -target  particle  filtering  method  that  provides  the 
baseline  to  understand  the  effectiveness  of  the  PHD  filter  to 
estimate  the  number  of  targets  and  localize  their  positions. 
Simulations  are  provided  in  Section  IV  to  characterize  the 
performance  of  the  PHD  filter  in  light  of  the  clairvoyant  track¬ 


ing  methods  and  ClusterTrack.  Two  sensor  model  matching 
cases  and  two  sensor  model  mismatching  cases  are  considered. 
Finally,  Section  V  provides  concluding  remarks. 

II.  Sensor  Model 

This  section  presents  two  sensing  models  that  describe  how 
the  interactions  of  the  targets  within  a  proximity  sensor  lead  to 
a  binary  measurement.  The  probabilistic  sensor  model  taken 
from  [2]  considers  how  the  signals  radiating  from  the  targets 
physically  manifest  into  a  measurement.  The  disc  model  is  an 
idealization  that  is  commonly  seen  in  the  literature.  While  we 
believe  that  the  probabilistic  model  is  a  better  representation 
of  how  measurements  emerge,  it  is  interesting  to  understand 
the  implication  of  both  models  on  tracking  performance  in  a 
proximity  sensor  network. 


A.  Probabilistic  sensor  model 

The  probabilistic  model  accounts  for  the  effects  of  sensor 
noise  so  that  a  sensor  will  report  (or  not  report)  a  detection 
based  on  some  probability  of  detection.  This  detection  proba¬ 
bility  goes  up  as  the  signal  from  the  targets  increases,  which 
happens  as  targets  move  closer  to  the  sensor.  The  probabilistic 
model  used  in  this  paper  is  taken  from  [2].  In  essence,  the 
model  assumes  that  target  signals  are  uncorrelated  so  that  the 
power  measured  at  the  sensor  is  the  sum  of  target  powers  at  the 
sensor  plus  measurement  noise.  Specifically,  the  instantaneous 
received  power  pi,  i.e.,  the  power  measurement  of  the  i-th 
sensor  at  a  given  point  in  time,  is  given  by 


where  po.fc  is  the  power  measured  at  a  reference  distance  ro,fc 
due  to  the  fc-th  target,  Nt  is  the  total  number  of  targets,  r\  k 
is  the  relative  distance  between  the  i-th  sensor  and  the  k- 
th  target,  a  is  the  attenuation  parameter  that  depends  on  the 
transmission  medium,  and  i>,;  is  additive  white  Gaussian  noise, 
i.e., 

Vi  ~  N  (pv,al)  . 


The  mean  and  variance  of  the  error  vt  is  derived  from  the 
zero  mean  measurement  noise  of  variance  a2  for  the  case 
that  the  measured  power  is  the  result  of  integrating  the  square 
of  the  measurements  over  L  samples.  As  shown  in  [2],  [22], 
integrating  over  L  samples  leads  to  pv  =  a2  and  a2  =  2 er4/L. 
This  paper  assumes  the  reference  distances  =  ro  and 
target  powers  pok  =  p0  for  all  K  targets,  and  r0  and  p0  are 
known.  As  discussed  later  in  this  subsection,  knowledge  of 
Po  enables  one  to  infer  the  number  of  targets  simply  from  the 
binary  detection  reports.  Clearly,  the  assumption  that  the  target 
signal  powers  are  all  equal  and  known  is  not  realistic,  and 
future  work  will  investigate  the  relaxation  of  this  assumption. 

The  i-th  sensor  measures  the  received  power  pi,  processes 
it  locally,  and  reports  a  single  binary  digit:  ‘  1  ’  for  the  presence 
of  one  or  more  targets  or  ’0’  for  the  absence  of  any  target. 
The  decision  follows  the  rule 


f  1  Pi  >  A, 
\  o  Pi  <  X. 


(2) 
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or  multiple  quiet  targets?  For  this  reason,  we  assume  in  this 
paper  that  poj,:  =  p0 .  Future  work  will  investigate  robustness 
of  this  assumption  when  the  targets  radiation  power  is  known 
to  fall  within  a  certain  range. 


B.  Disc  sensor  model 

Many  papers  exploiting  proximity  sensors  assume  the  disc 
sensor  model,  where  a  sensor  always  detects  one  or  more 
targets  if  they  are  within  the  detection  range  reff  and  never 
detects  any  target  if  all  are  further  than  ref  f  from  the  sensor: 


Fig.  1.  Pd  vs.  target  distance,  where  Pfa=  0.001,  a  =  2,  a  =  0.5,  po  =  P^is,  [Vi,  •  •  .  ,y/cl)  =  j  ^  ^  S^!  1  eff’  (5) 

3000,  ro  =  lm,  and  L  =  100.  1  0  otherwise. 


The  probability  of  false  alarm  Pfa  is  the  probability  that  z, 
exceeds  the  threshold  when  Nt  =  0.  For  a  desired  value  of 
Pfa,  the  threshold  A  is  computed  as: 

A  =  avQ~1(Pfa)  +  nv,  (3) 

where  Q(-)  is  the  Q- function 

Q(x)  =  _ _  /  e~f~  dt. 

V2t tJx 

Furthermore,  the  probability  of  detection  P,i  for  the  chosen 
threshold  (or  Pfa)  at  the  sensor  is  a  function  of  the  positions 
of  the  Nt  targets,  i.e.,  y/;:,  relative  to  the  location  of  the  sensor, 
i.e.,  s,  as  given  by 


The  disc  model  ignores  the  physical  process  as  discussed 
in  the  previous  subsection  that  enables  sensors  to  detect  the 
targets.  For  instance,  the  presence  of  multiple  targets  extends 
the  detection  range  of  the  sensor  due  to  the  superposition 
of  the  target  energy.  As  discussed,  it  is  the  superposition  of 
power  property  that  allows  one  to  estimate  the  number  of 
targets.  For  the  disc  model,  the  number  of  targets  is  always 
ambiguous.  Consider  a  proximity  sensor  network  with  one 
cluster  of  detections.  The  cluster  could  be  the  result  of  one 
target  or  hundreds  of  targets  within  the  same  vicinity.  Beyond 
the  fact  that  the  disc  model  is  an  idealized  representation  of  the 
behavior  of  the  sensor,  the  inability  of  the  disc  sensing  model 
to  disambiguate  the  number  of  targets  is  a  severe  weakness 
that  can  affect  the  performance  of  tracking  algorithms  based 
on  this  model  as  demonstrated  in  Section  IV. 


pd(s,  [yi,---,yivt]) 


Prob(zi  =  1 1  targets), 
Prob(pi  >  A | targets), 


rO,k 

rk 


(4) 


where  rp  =  ||yfc  —  s||.  Fig.  1  illustrates  Pd  as  a  function 
of  the  target  distance  for  various  numbers  of  targets.  For  the 
case  of  a  single  target,  the  figure  demonstrates  that  Pd  as  given 
by  (4)  behaves  as  a  sigmoid  function  that  steps  down  from  a 
value  near  one  to  a  value  near  Pfa  over  a  transition  region 
from  100m  to  200m. 

For  multiple  targets,  the  addition  of  the  target  powers  as 
seen  in  (4)  does  extend  the  sensing  range  of  the  sensors  as 
demonstrated  in  Fig.  1.  This  fact  actually  allows  for  one  to 
infer  the  number  of  targets  when  looking  at  the  distribution  of 
detection  reports  in  a  mesh  network  of  proximity  sensors.  If 
the  targets  are  far  apart,  one  would  expect  to  see  separate 
clusters  of  sensors  reporting  detections  around  each  target. 
When  the  targets  are  near  each  other,  only  one  cluster  emerges. 
The  size  of  the  cluster  increases  as  more  targets  are  in  the 
same  vicinity  due  to  the  addition  of  powers  in  (4),  and  it  is 
the  size  of  the  cluster  that  provides  the  clue  to  the  number  of 
targets.  The  ability  to  estimate  the  number  of  targets  depends 
on  having  a  good  handle  on  the  target  radiation  power  p0,fc-  If 
Po.fc  is  n°t  known,  then  the  number  of  targets  within  the  cluster 
of  sensing  nodes  is  ambiguous,  e.g.,  is  there  one  loud  target 


III.  Multiple  Target  Tracking 

The  multiple  target  tracking  approach  for  the  binary  proxim¬ 
ity  sensor  network  proposed  in  this  paper  exploits  the  theory  of 
finite  point  processes  (or  random  finite  set  statistics  (FISST)). 
Rigorous  derivations  of  multiple  target  tracking  from  FISST 
can  be  found  in  [11],  [13],  [18].  In  essence,  at  a  given  snapshot 
in  time,  the  set  of  all  targets  is  modeled  as  a  realization  of 
a  finite  point  process.  The  tracking  approach  simply  updates 
statistics  about  the  posterior  distribution  of  the  finite  point 
process  given  the  sensor  measurements. 


A.  Finite  Point  Process  Basics 

For  the  purposes  of  this  paper,  a  finite  point  process  is  a 
random  process  whose  realization  is  a  set  of  points  within  a 
d-dimensional  state  space  Kd.  The  realization  can  be  viewed 
as  a  sample  from  a  two  step  process  where  the  number  of 
objects  n  is  sampled  from  the  probability  mass  function  (pmfj 
fn(n),  and  then  the  n  d-dimensional  states  are  sampled  from 
the  joint  probability  density  function  (pdf)  fx i„(xi, . . .  ,xn|?t) 
that  is  conditioned  on  the  number  of  targets.  The  overall  state 
of  a  realization  of  a  finite  point  process  is  the  number  of 
points  and  the  individual  state  vectors  describing  each  point. 
Thus,  the  finite  point  process  state  is  a  variable  dimension  state 
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X  =  [  nx  xf  ...  x^  ]  of  size  ( nxd+  1)  x  1.  The  composed  of  a  prediction  and  a  likelihood  update  stage.  The 

probability  density  for  the  finite  point  process  state  is1  prediction  step  starts  with  the  PHD  Dt_i|t_i(x)  of  the  prior 

multiobject  density,  i.e.,  the  density  of  the  finite  point  process 
f(X)  =  /( nx ;x1;  ...,x„x)  =  fx  |n(Xl, . . . ,  xnx  \nx)fn(nx),  (6)  given  measurement  up  to  time  t-  1.  We  refer  to  Z>t_1,t_1(x) 

as  the  prior  PHD.  The  prediction  step  determines  the  PHD 
which  is  referred  to  as  the  multiobject  density  function  for  £)t |t_1  (x)  for  the  density  of  the  finite  point  process  as  it 

the  finite  point  process.  As  any  probability  density,  the  integral  evolves  to  time  t  We  refer  to  as  the  predicted  PHD. 

over  all  possible  realizations  is  one.  Given  that  the  number  of  The  evolution  considers  target  motion,  death,  and  birth  models 
arguments  within  the  density  differs  based  on  the  number  of  for  the  points  in  [11],  Mahler  proved  a  simple  relationship 
realized  points,  the  integration  of  all  realizations  is  between  the  prior  and  predicted  PHDs  (see  (75)  in  [11]).  This 

work  does  not  consider  the  target  death  and  birth  models,  and 

/T5L  r  JFL  without  these  models,  the  relationship  simplifies  to 

f(X)SX=  £  /  f{nx-,xi,...,xnx)dxi"-dxnx=  £  /"(«*)  =  !• 

nx=°  A|t-i(x)  =  J  /t|t-i(x|w).Dt_1|t_1(w)dw,  (9) 


f(X)  =  f(nx;x  i, . .  -,x„x)  =  fx  |n(xi, . .  .,xnx\  nx)f„(nx),  (6) 

which  is  referred  to  as  the  multiobject  density  function  for 
the  finite  point  process.  As  any  probability  density,  the  integral 
over  all  possible  realizations  is  one.  Given  that  the  number  of 
arguments  within  the  density  differs  based  on  the  number  of 
realized  points,  the  integration  of  all  realizations  is 


A  useful  statistics  to  summarize  the  multiobject  density  is 
the  probability  hypothesis  density  (PHD),  which  represents  the 
support  of  the  density  for  a  particular  target  state  over  all 
possible  target  count  hypotheses,  i.e., 


£>(x)  =  £  /  -x) )  /(nT;yi>---’y«y)dyi 


Therefore,  the  PHD  represents  the  density  of  the  expected 
number  of  objects  at  a  particular  target  state  x;  that  is,  the 
integral  of  the  PHD  in  the  target  space  represents  the  expected 
number  of  objects.  In  short,  the  PHD  filter  not  only  provides 
the  density  of  multiple  targets  but  also  the  estimates  of  the 
number  of  targets  in  the  scene. 

B.  Track  Filtering 

Single  target  track  filtering  basically  transforms  a  posterior 
target  density  for  time  t  —  1  into  a  predicted  (or  prior)  density 
for  time  t  through  the  use  of  a  stochastic  target  motion 
model.  Then  the  prior  density  is  transformed  to  a  posterior 
density  through  Bayes’  rule  in  the  likelihood  (or  measurement) 
update  stage.  Similarly,  it  is  possible  to  perform  multiple 
target  track  filtering  using  the  multiobject  density.  A  particle 
implementation  of  such  a  Bayes  multitarget  filter  is  presented 
in  [13].  The  main  problem  with  the  Bayes  multitarget  filter 
is  that  the  dimensionality  of  the  multiobject  density  function 
is  not  bounded,  and  as  the  number  of  targets  increases,  the 
support  of  the  multiobject  density  migrates  to  portions  of  the 
density  representing  increasingly  higher  dimensionality. 

The  probability  hypothesis  density  (PHD)  filter  was  intro¬ 
duced  by  in  [11]  as  a  means  to  control  the  curse  of  dimension¬ 
ality.  Instead  of  tracking  the  full  multiobject  distribution  over 
time  for  the  finite  point  process,  the  filter  simply  tracks  the 
PHD  as  defined  in  (8),  which  is  a  first  order  moment  of  the 
distribution.  Just  like  any  other  track  filter,  the  PHD  filter  is 

1  As  defined,  f(X)  is  virtually  the  same  as  —  Jn({xi, . . . ,  xn})  where 
jn  (-)  is  the  Janossy  density  [11],  For  the  sake  of  notational  convenience  in 
the  development  of  the  particle  filter,  the  symmetry  property  of  the  Janossy 
density  is  not  maintained  in  f(X),  e.g.,  the  density  (21)  is  technically  not 
symmetric.  Forcing  the  multiobject  density  to  be  symmetric  in  the  particle 
filter  development  is  straightforward  but  requires  additional  computational 
resources,  e.g.,  memory,  and  arrives  at  the  same  PHD  due  to  collapsing  the 
density  (see  (23))  and  the  symmetry  of  the  likelihood  given  by  (13). 


where  /t|i_1(x|w)  is  the  density  of  a  point’s  state  x  at  time  t 
conditioned  on  its  state  w  at  time  t  —  1.  The  specific  form  for 
/t|i_i(x|w)  represents  the  motion  model.  The  motion  model 
evokes  a  Markov  assumption  for  the  movement  of  points,  and 
it  assumes  targets  move  independently  of  each  other. 

The  likelihood  update  stage  evokes  Bayes’  rule  to  determine 
cfyhePHD  of  the  posterior  multiobject  density  given  the  mea¬ 
surements  up  to  time  t.  In  general,  the  resulting  PHD  Dt\t(x), 
which  we  refer  to  as  the  posterior  PHD,  cannot  be  uniquely 
determined  solely  by  the  predicted  PHD.  In  other  words,  two 
mutiobject  density  functions  that  have  the  same  PHD  will  not 
in  general  have  the  same  PHD  after  the  likelihood  update.  To 
ground  the  density,  the  PHD  filter  approximates  the  predicted 
multiobject  density  with  the  density  of  a  Poisson  point  process 
(PPP).  As  show  in  Theorem  4  in  [11],  the  best  Poisson 
approximation  in  the  sense  of  the  Kullback-Leibler  divergence 
is  the  Poisson  point  process  whose  intensity  function  is  the 
predicted  PHD,  i.e.. 


/tit-iW  =  .ft\t-i(n-,x1, . . . . 


0=^r 

n! 


riAn-iM 


Given  the  density  of  the  measurements  conditioned  on  the 
realization  of  the  point  process,  Bayes’  rule  leads  to  the 
posterior  multiobject  density 


/t|t(n;  x1; . . . ,  xn)  =  ^fz\x(z\X)^ 


II  Al, -!(**) 


where  C  is  the  constant  so  that  the  posterior  integrates  to  one. 


c=£  nAit-i(xi) 


Then,  the  posterior  PHD  Dt|t(x)  is  calculated  by  applying 
(8)  over  (11).  Finally,  Dt\t(x)  is  treated  as  the  PHD  for  the 
prior  multiobject  density  for  the  next  round  of  prediction  and 
measurement  updates. 

In  the  traditional  use  of  the  PHD  as  proposed  by  Mahler,  the 
measurements  are  modeled  as  a  point  process  whose  points  are 
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a  mixture  of  a  clutter  PPP  and  points  spawned  by  targets.  In 
this  case,  fz\x  is  a  multiobject  density  function  conditioned  on 
the  realization  of  the  target  process.  Furthermore,  an  analytical 
relationship  exists  between  the  posterior  and  predicted  PHDs 
so  that  the  use  of  Bayes’  rule  over  the  PPP  approximation  for 
the  multiobject  density  is  implicit.2  As  a  result,  the  particle 
implementation  of  the  traditional  PHD  filter  only  requires 
particles  to  represent  the  PHD.  In  general,  we  view  the  PHD 
filter  as  the  process  of  performing  the  prediction  and  likelihood 
updates  for  the  tracking  of  the  PHD  of  the  finite  point  process 
conditioned  on  the  observations  via  (9),  (11),  and  (8).  When 
the  measurements  are  a  mixture  of  a  clutter  PPP  and  target 
spawned  measurements,  then  the  PHD  filter  is  the  familiar 
implementation. 

For  the  application  of  multiple  target  targeting  using  binary 
sensor  measurements  over  a  2-D  field,  the  target  state  is  a  4-1? 
vector  representing  the  2-D  position  y  and  2-D  velocity  v  of 
the  target,  i.e.,  x  =  [  yT  vT  ]  and  d  =  4.  The  discretized 
continuous-time  white  noise  acceleration  (DCWNA)  model 
is  employed  for  the  motion  model  [23].  Each  binary  sensor 
measurement  has  a  pmf  given  by  (4)  or  (5)  for  the  probabilistic 
or  disc  model,  respectively.  The  measurement  model  also 
assumes  the  sensor  measurements  are  statistically  independent. 
For  a  proximity  sensor  network  of  Ns  sensors,  the  likelihood 
fz\x  f°r  measurements  z  £  {0,l}‘Ys  is  given  by 

Ns 

fz\x(*\x)  =  n  Pd(sm,yym  (i  -  Pd^y))1-^ , 

m=  1 

(13) 

where  y  =  [yi, _ ,  y„]  is  extracted  from  X,  is  the  location 

of  the  2-th  sensor,  and  z  is  the  collection  of  binary  measure¬ 
ments  at  the  current  time.  Note  that  the  likelihood  is  a  pmf 
evaluated  at  the  values  of  the  reported  sensor  measurements. 
From  (4)  (or  (5))  and  (13),  it  is  clear  that  the  likelihood  is 
symmetric  with  respect  to  any  permutation  of  the  target  state 
indices  within  X.  To  the  best  of  our  knowledge,  the  calculation 
of  the  posterior  PHD  from  the  posterior  multiobject  density 
after  the  likelihood  update  does  not  simplify  as  a  simple 
analytical  transformation  of  the  prior  PHD. 


C.  Particle  PHD  Filter  for  Binary  Sensors 

Unlike  the  particle  implementation  of  the  traditional  PHD 
filter,  the  particle  implementation  of  the  general  PHD  filter  re¬ 
quires  particles  to  sample  the  PHD,  and  other  particle  derived 
from  the  PHD  particles  to  sample  the  predicted  multiobject 
density.  The  basic  steps  of  the  proposed  particle  PHD  filter 
are  to  1)  propagate  the  PHD  particles  via  prediction  update  in 
(9),  2)  sample  the  predicted  multiobject  density  from  the  PPP 
approximation  in  (10)  to  form  multiobject  density  particles, 

3)  compute  the  likelihood  update  in  (11)  using  the  multiobject 
density  particles  as  an  approximation  of  the  actual  density,  and 

4)  collapse  the  multiobject  density  particles  into  PHD  particles 
through  (8).  This  section  explains  how  to  perform  these  steps. 

-[11]  uses  probability  generating  functionals  to  derive  the  PHD  likelihood 
update  equation  where  the  usage  of  Bayes’  rule  via  (11)  is  also  implicit. 


To  begin,  the  prior  PHD  at  time  i  —  1  is  approximated  by 
P  single  target  particles  x^1}^=1  so  that 

p 

P>t- i|t-i(x)  «  w^S (x  -  x^),  (14) 

p=i 

where  S(-)  is  the  Dirac  delta  function.  Then,  the  estimated 
number  of  targets  is  the  integral  of  the  PHD  so  that 

iVt_t  «  5>t-f  (15> 

p=i 

To  compute  the  predicted  PHD,  the  p-th  particle  is  sampled 
from  the  proposal  density  x*^  ~  /t|t_1(x|x^1).  Given  the 
DCWNA  motion  model,  the  predicted  particles  are  obtained 
by  diffusing  the  prior  particle  such  that 

x*(p)  =  Fx|!_)1  +  Q5V(p),  (16) 

where  vfp)  is  sampled  from  a  white  Gaussian  random  number 
generator,  and  F  and  Q  are  the  state  transition  and  process 
noise  covariance  matrices,  respectively  (see  [23]).  In  other 
words,  the  proposal  density  is  Gaussian  with  mean  Fx^\  and 
covariance  Q.  Since  the  proposal  density  matches  the  actual 
motion  model,  the  predicted  PHD  is  simply 

A|t-i(x)  «X!wt^(x“Xt(p)).  (17) 

p=i 

Note  that  since  the  target  prediction  step  does  not  incorporate  a 
target  birth  and  death  model,  the  predicted  PHD  still  integrates 
to 

To  perform  the  likelihood  update,  samples  from  the  PPP 
approximation  of  the  predicted  multiobject  density  given  by 
(10)  must  be  generated  where  the  predicted  PHD  is  (17). 
The  two  stage  PPP  sampling  approach  as  explained  in  [24] 
is  employed.  For  the  p- th  multiobject  particle,  the  number  of 
hypothesized  targets  is  sampled  from  a  Poisson  distribution 
with  mean  Nt_ i,  that  is,  n ^  ~  Poisson{Nt_\).  Then,  each 
of  the  states  of  the  targets  is  sampled  from  the  normalized 
predicted  PHD  given  by  (17)  divided  by  Nt- 1-3  This  is 
equivalent  to  selecting  each  of  the  target  states  to  be  one  of  the 
P  PHD  particles  where  the  particle  to  select  is  sampled  from 

w&> 

the  multinomial  distribution  given  by  J  1  for  p  =  1 , ,P. 
The  proposed  method  uses  this  PPP  sampling  approach  to 
generate  Pnex  multiobject  particles  where  nex  >  1  is  a  user 
supplied  parameter. 

When  a  modest  number  of  particles  are  used  to  approximate 
the  predicted  PHD,  it  is  possible  for  the  particles  to  drift 
away  from  the  true  target  locations  due  to  an  imperfect  motion 
model.  To  regularize  the  process,  it  is  helpful  to  generate  some 
innovative  multiobject  particles  by  sampling  from  likely  target 
locations  given  the  current  set  of  measurements.  For  these 
particles,  the  number  of  targets  represented  by  a  single  multi¬ 
object  particle  is  still  sampled  from  a  Poisson  distribution  with 
mean  Nt-i-  The  target  states  are  generated  from  a  proposal 
density  q(- |z)  that  is  conditioned  on  the  current  measurements. 

^Normalized  to  integrate  to  one. 
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This  proposal  density  could  be  the  scaled  likelihood  using 
the  Metropolis-Hastings  method  [25],  but  such  a  method  is 
of  high  computational  complexity.  In  this  paper,  the  proposal 
density  for  the  position  components  of  the  target  states  is  a 
rectangular  bounding  box  that  encompasses  all  sensors  that 
report  detections.  The  velocity  components  of  the  target  states 
are  set  to  zero.  The  bounding  box  represents  a  computationally 
efficient  proposal  density  that  approximates  the  more  costly 
method  of  incorporating  the  feasible  target  area  (FTA)  used  in 
[20],  [21]. 

Overall,  J  innovative  multiobject  particles  are  generated, 
and  Pex  =  Pnex  +  J  multitarget  particles  represent  the 
predicted  target  density 


ft\t- i{X) 


E  *(*;*“). 

p= i 


(18) 


where 

~S(X-,Xt{p)) 


0  ^  if  \x\^\x^\ 

nr=i  <Kx-xm)  otherwise 


(v) 

is  the  equivalent  of  a  delta  function  centered  at  Xt  over  the 
multiobject  space,  and  rif^  and  x^p}  are  extracted  from  X ^ .  It 
is  not  difficult  to  verify  that  <5(-)  satisfies  the  sampling  property 

j  f(X)5(X;  Xt{p))SX  =  f(Xt{p))  (19) 


so  that 

f(X)~S(X-  Xt(p))  =  f(X^)6(X-  Xtip)).  (20) 


Also,  the  weights  are  not  explicit  in  (18)  since  they  are 
equal  to  one  due  to  the  fact  the  particles  are  sampled  from 
the  desired  density.  In  effect,  the  generation  process  creates 
particles  sampled  from  a  PPP  with  an  intensity  (or  PHD)  of 

D(x)  =  +  ^-A-(x), 

*  ex  * ex 

where  D,.(x)  represents  the  portion  of  the  regulated  PHD  due 
to  the  innovative  particles.  It  is  uniform  over  the  bounding  box 
surrounding  the  sensors  reporting  detections,  and  it  integrates 
out  to  Nt_ i  so  that  the  estimated  number  of  targets  due  to 
regularization  does  not  change.  The  ratio  J  over  Pex  represents 
the  importance  the  regularization  places  over  on  current  mea¬ 
surements  relative  to  the  prior  measurements.  When  J  =  0, 
the  predicted  distribution  is  the  normal  predicted  distribution. 
On  the  other  hand,  when  J  =  Pex  the  predicted  distribution 
completely  ignores  the  prior  measurements.  A  moderate  value 
of  J  prevents  the  set  of  particles  to  lose  track  of  the  target  due 
1)  to  either  too  few  particles  to  properly  represent  the  density 
or  2)  to  an  overly  simplistic  dynamical  model. 

The  likelihood  update  passes  the  density  given  by  (18) 
through  (11)  and  (12)  using  the  likelihood  given  by  (13).  As  a 
result  of  the  sampling  property  (19)-(20),  the  posterior  density 
is 

pex 

ftlt(X)^Y,wfp'S(X’XtP)),  (2D 

P=  1 


where  the  weights 

Wp  = -^fz\x{zt\Xtip)),  (22) 

zt  are  the  binary  measurements  collected  at  time  t,  and  C  is 
a  constant  so  that  the  weights  sum  to  one. 

The  final  step  of  the  PHD  filter  is  to  collapse  the  posterior 
density  back  down  to  the  posterior  PHD  by  inserting  (21)  into 
(8).  As  a  result, 

Pex  n(p) 

A|t(x)  «  E^E^x-xg),  (23) 

p—  1  i—1 

~  Ew*(p)<5(x_x*(p))’ 

p= i 

where  P*  =  Y^p=i  n<J>\  and  there  is  a  one-to-one  corre¬ 
spondence  between  the  PHD  particles  x( !/,;  and  multiobject 
particle  components  xj  (  .  Each  multiobject  particle  component 
leads  to  a  PHD  particle,  and  when  the  fc-th  PHD  particle 
originated  from  the  A-th  component  of  the  m-th  multiobject 
particle,  then  x^fc-*  =  xj1"^  and  uE'  =  w^.  The  expression 
for  the  conversion  of  multiobject  particles  to  PHD  particles 
has  appeared  in  [13]  (see  Section  2E)  for  generating  target 
state  estimates  in  implementation  of  the  Bayes  multitarget 
filter.  Integration  of  (23)  leads  to  a  target  count  estimate  of 
Nt  =  E3  r^  wl  =  Finally,  the  PHD  is 

resampled  down  to  P  particles  so  that  the  number  of  samples 
does  not  continue  to  grow  after  each  likelihood  update.  To  this 
end,  the  P  particles  are  selected  from  the  Pex  particle  using 
sampling  with  replacement  from  the  multinomial  distribution 

w *(p) 

—A — .  After  this  final  resampling  stage, 

A|tW  «  E  w[p)6{x  -  xiP))>  (24) 

P=  1 

where  w[p '  =  Nt/P.  For  the  next  stage  of  prediction  and 
measurement  updates,  (24)  replaces  (14)  as  the  prior  PHD. 

The  implementation  of  the  PHD  particle  filter  at  any  time  is 
summarized  in  Fig.  2:  Start  with  the  prior  PHD  particles  (Step 
1),  form  the  predicted  PHD  particles  via  the  target  motion 
model  (Step  2),  expand  to  the  particle  representation  of  the 
PPP  multiobject  density  (Steps  3-5),  perform  the  measurement 
(or  likelihood)  update  (Step  6),  and  calculate  the  posterior 
PHD  (Step  7-10).  After  the  P  particles  for  Dt |t (x)  are 
obtained,  the  A- means  method  is  used  to  group  the  particles 
into  round(-/Vt)  clusters.  The  centroids  of  these  clusters  are 
used  as  the  estimated  target  states. 

The  particle  implementation  of  the  PHD  filter  in  Fig.  2  is 
able  to  estimate  the  number  of  targets  without  any  target  birth 
and  death  model  due  to  the  expansion  of  the  PHD  particles  to 
multiobject  particles.  Specifically,  Step  3  means  that  the  PHD 
does  explore  hypotheses  of  various  numbers  of  targets.  Then, 
the  likelihood  update  is  able  to  emphasize  hypotheses  whose 
target  counts  better  explain  the  measurements  z  at  the  given 
snapshot  in  time.  It  is  these  hypotheses  that  mostly  determine 
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1)  Given  P  prior  PHD  particles  where  Nt-i  =  Sp=i  wt- 1- 

2)  Obtain  P  predicted  PHD  particles:  x.^p'>  =  +  Q^v^p\p  =  1, 2, . . . ,  P,  where  is  a 

zero  mean  white  Gaussian  noise  vector. 

3)  Obtain  Pex  Poisson  samples  with  mean  Nt-i ■  n ^  ~  Poisson(Nt_i),  k  =  1, 2, ... ,  Pex ,  where 

Pex  J  • 

4)  Obtain  Pnex  predicted  multiobject  particles  {xj^, . . . ,  x|fc^fc) }£2jx  where  the  fc-th  multiobject 


„(p) 


particle  consists  of  n target  states  drawn  (with  replacement)  from  {x^ p  }£L1  with  probability 


(p) 

Ut- 1 


N,. 


5)  Obtain  J  innovative  multiobject  particles  {xj^, . . . ,  \k=Pn  J+  t  where  the  fc-th  multiob¬ 

ject  particle  consists  of  n<k)  target  state  samples  drawn  uniformly  from  a  rectangular  box  that 
bounds  all  the  detecting  sensors. 

6)  Obtain  the  weights  {wl} {ffi  for  the  Pex  multiobject  particles  {x^ 
and  then  normalize  the  weights  to  sum  to  one. 

7)  Convert  Prx  multiobject  particles  to  P*  PHD  particles  at  time  t  by  (21)  to  form 


,(*) 

H.i  >  ■ 


{(fc) 

t,n 


by  (13), 


(  *(k)  *trej-|P' 

{«V  }fe=l-  , 

8)  Estimate  the  number  of  objects  at  time  t:  Nt  =  Y^k= l  n^wl- 

9)  Obtain  P  posterior  PHD  particles  drawn  (with  replacement)  from  {xj(p^}^=1  with  probability 


ex 

p* 


JV. 


10)  Weigh  each  PHD  particle  x^  uniformly  as  w\p>  =  Nt/P. 

11)  Let  t  :=  t  +  1  and  goto  Step  #1. 


,(p) 


Fig.  2.  Steps  of  the  particle  implementation  of  the  PHD  filter  for  proximity  sensor  networks. 


the  posterior  PHD  particles.  Nevertheless,  the  particle  imple¬ 
mentation  could  be  enhanced  by  explicit  target  birth  and  death 
models.  This  is  left  for  future  investigation. 


D.  Initialization 

The  PHD  filter  begins  with  an  initial  set  of  particles  at  time 
t  =  0  to  approximate  Z?0|0(x).  For  most  of  the  simulations  in 
Section  IV,  we  generated  P  particles  where  the  location  states 
are  uniformly  sampled  over  the  surveillance  region  covered 
by  the  sensors,  and  the  velocity  states  are  sampled  from  a 
Gaussian  with  mean  zero  and  a  =  10.  The  weight  associated 
to  each  particle  is  initially  set  to  N/P  where  N  is  a  prior 
estimate  of  the  number  of  targets.  In  the  simulations,  N  is 
set  up  to  be  0.5  for  the  PHD  filter.  4  For  the  clairvoyant 
particle  filters  in  Section  III-E  and  Section  III-F,  N  is  the  true 
number  of  targets  because  we  assume  the  number  of  targets  is 
known  a  priori.  To  determine  the  sensitivity  of  the  PHD  to  the 
initialization,  some  of  the  simulations  generated  the  location 
states  by  restricting  these  values  to  be  drawn  uniformly  within 
the  rectangular  bounding  box  for  sensors  reporting  detections. 

E.  Clairvoyant  multitarget  particle  filter 

Unlike  the  PHD  filter,  the  particle  implementation  of  mul¬ 
titarget  Bayes  filter  requires  explicit  target  birth  and  death 
models.  For  example,  when  the  scene  contains  N  targets, 
the  multiobject  particles  that  survive  will  correspond  to  N 
targets.  Without  explicit  birth  and  death  models,  when  more 
targets  emerge  or  if  targets  leave  the  scene,  there  is  no 

4This  is  to  show  that  despite  the  mismatch  in  the  initialization,  the  target 
count  can  quickly  converge  to  the  true  value. 


mechanism  for  multiobject  particles  corresponding  to  target 
counts  different  than  N  to  emerge.  To  perform  a  comparison 
of  the  PHD  filter  with  the  multitarget  particle  filter  in  this 
paper  where  birth/death  models  are  not  considered,  a  clair¬ 
voyant  multiobject  particle  filter  is  used.  This  allows  one  to 
understand  the  effectiveness  of  the  PHD  filter  to  track  multiple 
targets.  The  particle  filter  is  clairvoyant  in  the  sense  that 
the  number  of  targets  N  is  known  a  priori.  Therefore,  no 
target  birth/death  model  is  needed.  The  performance  of  the 
filter  indicates  the  effectiveness  of  standard  processing  when 
the  number  of  targets  need  not  be  estimated.  Effectively,  the 
A? -target  clairvoyant  particle  filter  (CPF)  is  a  particle  filter 
implementation  of  the  multitarget  sequential  Bayesian  filter 


when  the  multiobject  density  function  f(nx ;  xy, . . . 


, )  has 


nonzero  support  only  for  nx  =  N.  Thus,  the  p- th  multiobject 


x(p)  x(p) 

JVj.  i  ...  -A., 


particle  is  X^p)  =  N  ^tl  ...  N 

Initially  at  time  t  =  0,  P  multiobject  particles  are  generated. 
These  particles  are  composed  of  N  4-state  subparticles  that 
collectively  form  a  sampling  of  fo\o(N",  xi, . . . ,  xjv).  For 
cases  where  targets  appear  or  disappear  at  certain  times,  the 
multiobject  density  function  is  adjusted  accordingly.  Given  P 
multiobject  particles,  the  evolution  of  the  particles  }^=1 

at  time  t  —  1  proceeds  as  follows: 


(r>) 

1)  Obtain  P  predicted  multiobject  particles:  ’  = 


r  (pUat 
lxt,i  Ji=l 
1,2, 


where  x 


2) 


3) 


P. 

Obtain  J  innovative 

r(p)\N 


Ip) 

t,i 


=  f'x 


(p) 

t—l,i 


+  = 


(v) 

multiobject  particles  Xt  = 


{xjv  })!]  where  the  p- th  multiobject  particle  consists 
of  N  target  states  drawn  uniformly  from  a  rectangular 
box  that  bounds  all  the  detecting  sensors. 

Obtain  the  weights  {tUp  }p=\  for  the  P  +  J  multiobject 
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particles  {xjp>  where  the  weight  for  the  p-th 
multiobject  particle  wjj  is  computed  by  (13). 

4)  Resample  the  P  +  J  multiobject  particles  to  obtain 
P  particles  {X^}^=1  where  X^  is  selected  with 

probability  -ff. 

5)  Let  t:=t+l  and  goto  Step  #1. 

The  components  (target  state  vectors)  in  the  multiobject 
particles  may  not  align  between  particles.  In  other  words, 
the  *-th  target  component  in  the  particle  xfm>  and  the  i-th 
target  component  in  the  particle  x[n’  when  m  ^  n  may 
not  correspond  to  the  same  target.  To  calculate  a  minimum 
mean  squared  error  (MMSE)  estimate  of  the  target  states,  the 
subparticles  in  each  particle  are  assigned  ground  truth  targets 
by  solving  the  2 -D  assignment  problem  [26].  That  is  each 
subparticle  is  matched  to  one  and  only  one  ground  truth  target 
and  vice  versa.  Finally,  the  weighted  average  of  all  subparticles 
associated  to  a  particular  ground  truth  target  is  computed  as  the 
state  estimate  for  that  target.  Clearly,  the  CPF  is  incorporating 
information  that  is  unknown  to  the  filter  to  provide  the  state 
estimates.  The  purpose  of  the  CPF  is  to  simply  baseline  the 
performance  of  the  PHD  method. 

F.  Clairvoyant  PHD 

In  the  CPF,  the  grouping  of  subparticles  into  particles 
affects  tracking  performance  because  this  grouping  is  always 
maintained  in  the  filter.  On  the  other  hand,  the  PHD  continues 
to  mix  up  target  state  particles  when  it  forms  the  multiobject 
density  from  the  PHD.  This  process  allows  for  innovative 
combinations  of  single-target  particles.  In  some  sense,  the 
PHD  filter  performs  a  crossover  step  similar  to  genetic  pro¬ 
gramming.  However,  the  PHD  is  also  estimating  N,  which  the 
CPF  does  not  have  to  do.  To  determine  the  effectiveness  of 
these  crossover  operations  without  the  burden  of  estimating 
N,  we  also  consider  a  modified  PHD  where  N  is  known. 
Specifically,  Step  #3  of  the  PHD  filter  in  Fig.  2  is  changed  so 
that  n =  N  for  k  =  1, ... ,  Pex.  We  refer  to  this  PHD  as 
the  clairvoyant  PHD. 

G.  Performance  Measures 

The  Wasserstein  [27]  and  optimal  subpattern  assignment 
(OSPA)  [28]  metrics  are  common  measures  of  performance 
for  PHD  track  filters.  They  combine  the  effects  of  position 
and  cardinality  errors  in  a  single  quantity.  However,  we  wish 
to  deconflict  these  two  types  of  errors.  The  cardinality  error 
is  simply  \Nt  —  Nt\  where  Nt  and  Nt  are  the  estimated  and 
actual  number  of  targets  in  the  scene,  respectively. 

For  the  target  state  errors,  we  modify  the  OPSA  metric  [28] 
to  consider  only  the  position  errors  and  not  the  cardinality 
errors.  This  performance  measure  operates  on  the  set  of 
estimated  target  states  X  and  the  set  of  ground  truth  target 
states  Y  representing  the  position  components.  We  refer  to 
the  modified  performance  measure  as  the  cardinality  agnostic 
OPSA  (CA-OPSA).  The  CA-OPSA  is  not  a  metric  in  the  strict 
sense  as  a  measure  of  zero  does  not  imply  X  =  Y  due  to 
possible  differences  in  the  sizes  of  sets  X  and  Y .  Formally, 


1  1*1  m 

min  —  EE  ^i,j  d(xi  iDj)-,  (25) 

*= 1  3= 1 

E,-,  Oi  <  1.  ES^<1, 

1*1  p'l 

E  E  = 

i= 1  0= 1 

S  =  min(|X|,  |Y|),  and  Cij  G  {0, 1}, 

where  d(xi,yj)  is  a  user  selected  distance  measure  between 
the  i- th  and  / -th  position  states  for  sets  X  and  Y,  respectively. 
In  this  paper,  d( ■,■)  is  the  standard  f'1  norm.  In  essence,  the 
CA-OPSA  represents  the  average  root  mean  squared  (RMS) 
error  for  the  2 -D  assignment  of  the  estimated  and  ground  truth 
target  states  such  that  one  and  only  one  estimate  is  assigned 
to  a  ground  truth  element  and  vice  versa.  If  |X|  >  in 
the  extra  target  estimates  are  not  assigned  to  ground  truth 
elements.  Likewise,  if  \X\  <  |Y|,  the  extra  ground  truth 
elements  are  not  assigned  to  target  estimates.  Note  that  other 
metrics  do  exist  that  are  relatively  insensitive  to  differences 
in  the  target  cardinality,  e.g.,  the  Hausdorff  distance  [27],  For 
the  sake  of  simplicity,  this  paper  only  considers  the  CA-OSPA 
to  characterize  the  localization  performance  of  the  proposed 
PHD  filter. 

IV.  Simulations 
A.  Mode  Matching/Mismatching  Cases 

The  first  set  of  simulations  considers  four  scenarios  where 
targets  move  through  a  surveillance  region  of  size  1km  x 
lkm  for  100  or  300  seconds.  Filter  updates  occur  in  snapshot 
intervals  of  one  second.  The  geometry  of  Scenario  A  is 
illustrated  in  Fig.  3.  In  this  scenario,  two  targets  travel  at 
a  constant  velocity  with  a  speed  of  about  lOm/s.  The  two 
targets  cross  paths  at  t  =  50s.  Fig.  4  provides  the  geometry 
for  Scenario  B.  In  this  scenario,  two  targets  are  maneuvering. 
The  two  targets  never  cross  paths,  but  they  do  converge  and 
diverge  with  a  closest  approach  occurring  around  t  -  37s. 
Fig.  5  illustrates  Scenario  C  where  three  targets  move  in  and 
out  of  the  surveillance  region  during  the  100  seconds.  Finally, 
Fig.  6  illustrates  Scenario  D  where  five  targets  move  in  and 
out  of  the  surveillance  region  during  the  300  seconds.  The 
lower  number  at  an  end  of  a  track  represents  the  starting  time 
of  the  track.  Likewise,  the  larger  number  represents  the  end 
time  of  the  track.  For  example,  the  first  track  starts  at  t  =  1 
and  ends  at  t  =  300,  the  second  starts  at  t  =  31  and  ends  at 
t  =  270,  and  the  third  starts  at  t  =  61  and  ends  at  t  =  210, 
etc.  As  a  result,  the  number  of  targets  varies  from  1  to  5  at 
any  given  time. 

We  believe  that  the  probabilistic  sensing  model  provides 
a  more  realistic  representation  than  the  disc  model  for  the 
behavior  of  a  proximity  sensor  to  its  environment.  However, 
we  are  interested  in  understanding  the  robustness  of  multi¬ 
target  tracking  when  considering  either  model.  Therefore,  the 
simulations  consider  both  models  for  actual  sensing  behavior 
and  for  insertion  into  the  tracker.  Overall,  there  are  four  cases. 
The  first  two  are  matching  cases:  1)  the  actual  sensing  model 


the  CA-OPSA  is 

D(X,Y)  = 

s.t. 
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Fig.  3.  Scenario  A:  (a)  Geometry  of  a  100-sensor  network  with  two  targets 
traveling  at  constant  velocities,  and  (b)  separation  distance  between  the  two 
targets  vs.  time. 


is  the  probabilistic  model  and  the  filter  assumes  the  same 
probabilistic  model,  and  2)  the  actual  sensing  model  is  the 
disc  model  and  the  filter  assumes  the  same  disc  model.  The 
final  two  are  mismatching  cases:  3)  the  actual  sensing  model  is 
the  probabilistic  model  but  the  filter  assumes  the  disc  model, 
and  4)  the  actual  sensing  model  is  the  disc  model  but  the  filter 
assumes  the  probabilistic  model. 

For  the  probabilistic  sensing  model,  the  threshold  is  set  so 
that  Pfa  =  0.001.  In  the  disc  model,  re//  =  117m,  which  is 
the  distance  that  makes  Pd  =  0.5  for  Pfa  =  0.001  in  the  prob¬ 
abilistic  model.  The  filter  parameters  are  set  to  be  P  =  1000, 
nex  =  1,  and  J  =  100  so  that  Pex  =  1100.  For  each  of 
the  four  scenarios,  we  generated  10  random  configurations 
of  100  proximity  sensors,  and  for  each  configuration,  we  ran 
10  Monte  Carlo  realizations  of  sensor  reports.  Overall,  each 
scenario  consists  of  100  simulations.  The  reported  cardinality 
and  RMS  position  errors  (see  (25))  for  a  given  snapshot  time 
are  averaged  over  the  100  simulations.  For  these  simulations, 
we  compared  the  proposed  PHD  method  against  the  two 
clairvoyant  filters.  To  the  best  of  our  knowledge,  no  other 
baseline  method  exists  for  proximity  sensors  that  can  employ 
both  the  probabilistic  and  disc  sending  models. 

Fig.  7  plots  the  average  RMS  position  errors  versus  time 
for  the  PHD  and  two  clairvoyant  filters  for  Scenario  A  over 
the  four  sensing  model  cases.  Likewise,  Figs.  8-  10  provide 
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Fig.  4.  Scenario  B:  (a)  Geometry  of  a  100-sensor  network  with  two  targets 
maneuvering,  and  (b)  separation  distance  between  the  two  targets  vs.  time. 


Fig.  5.  Scenario  C:  Geometry  of  a  100-sensor  network  with  three  targets 
that  move  in  and  out  of  the  scene  at  different  times. 

similar  plots  for  Scenarios  B,  C  and  D,  respectively.  The 
relative  performances  of  the  PHD  compared  to  clairvoyant 
filters  are  consistent  over  the  four  scenarios.  Regardless  of 
the  actual  sensing  model,  the  trackers  that  incorporate  the 
probabilistic  model  into  the  measurement  likelihood  update 
provide  better  localization  accuracy  than  those  that  use  the 
disc  model.  The  probabilistic  model  accommodates  errors  in 
the  multiobject  particle  better  than  the  disc  model  because 
perturbations  in  the  particles  for  the  probabilistic  model  leads 
to  a  gradual  degradation  of  the  likelihood;  whereas,  perturba- 
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Fig.  6.  Scenario  D:  Geometry  of  a  100- sensor  network  with  five  targets  that 
move  in  and  out  of  the  scene  at  different  times. 


tions  in  the  particles  can  lead  to  sudden  one  to  zero  jumps 
in  the  likelihood  for  the  disc  model.  These  properties  might 
explain  the  better  performance  by  the  probabilistic  model. 
When  the  filters  incorporate  the  probabilistic  models  (Cases  1 
and  4),  the  filters  perform  similarly  for  Scenarios  A  and  B.  It 
takes  about  10  seconds  for  the  filters  to  reach  a  steady  state 
after  initialization,  and  the  localization  performance  degrades 
only  slightly  when  the  target  separation  is  small.  Overall,  the 
localization  performances  of  the  filters  when  implementing  the 
probabilistic  model  are  robust  when  the  actual  sensor  behavior 
follows  the  disc  model. 

For  the  matching  probabilistic  case,  the  clairvoyant  PF  is 
clearly  the  best,  and  the  clairvoyant  PHD  is  slightly  better  than 
the  PHD.  For  Scenario  D,  the  clairvoyant  PHD  is  better  than 
the  PHD  when  the  number  of  targets  is  small,  but  comparable 
when  the  number  of  targets  is  large  when  tracking  is  steady. 
When  the  filters  incorporate  the  disc  model,  the  performance 
of  the  CPF  is  much  poorer;  especially  for  the  mismatch  case 
(Case  3).  For  Case  2,  the  CPF  performance  is  good  when 
the  targets  are  near  each  other,  and  the  performance  degrades 
significantly  when  targets  are  far  apart.  This  is  an  artifact  of 
both  the  proposal  function  for  innovative  particles  and  the  lack 
of  the  crossover  effect  that  occurs  in  the  PHD  implementation. 
It  is  possible  that  the  weights  for  all  particles  become  zero  as  to 
be  discussed  later  in  Section  IV-D.  When  this  happens,  all  zero 
weighted  particles  are  then  given  equal  positive  weights.  As  a 
result,  this  weighting  process  then  cannot  distinguish  between 
particles  whose  subparticles  are  closer  to  the  detecting  sensors 
than  others.  This  problem  is  exacerbated  when  the  sensors 
actually  follow  the  probabilistic  model  where  false  alarms  and 
missed  detections  can  occur.  On  the  other  hand,  the  crossover 
effect  in  the  PHD  makes  the  probability  much  smaller  that  all 
particles  have  zero  weights.  This  might  explain  the  improved 
performance  for  the  PHD  methods. 

For  Scenario  C  and  D,  all  filters  exhibit  higher  errors 
when  new  target  appear  and  disappear  in  the  scene,  but  the 
performances  of  these  filters  improve  when  the  targets  are 
simply  moving.  Note  the  CPF  reinitializes  when  the  new 
targets  enter  or  leave  the  scene  while  the  PHD  filters  still 
propagate  the  previous  particles.  The  estimation  error  when 
the  number  of  targets  is  large,  is  higher  than  that  when 


Fig.  7.  Scenario  A:  CA-OSPA  errors  for  (a)  matching  probabilistic  models, 
(b)  matching  disc  models,  (c)  probabilistic  sensing  model  with  filter  assuming 
the  disc  model,  and  (d)  disc  sensing  model  with  filter  assuming  the  proba¬ 
bilistic  model. 
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Fig.  8.  Scenario  B:  CA-OSPA  errors  for  (a)  matching  probabilistic  models, 
(b)  matching  disc  models,  (c)  probabilistic  sensing  model  with  filter  assuming 
the  disc  model,  and  (d)  disc  sensing  model  with  filter  assuming  the  proba¬ 
bilistic  model. 
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Fig.  9.  Scenario  C:  CA-OSPA  errors  for  (a)  matching  probabilistic  models, 
(b)  matching  disc  models,  (c)  probabilistic  sensing  model  with  filter  assuming 
the  disc  model,  and  (d)  disc  sensing  model  with  filter  assuming  the  proba¬ 
bilistic  model. 
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Fig.  10.  Scenario  D:  CA-OSPA  errors  for  (a)  matching  probabilistic 
models,  (b)  matching  disc  models,  (c)  probabilistic  sensing  model  with  filter 
assuming  the  disc  model,  and  (d)  disc  sensing  model  with  filter  assuming  the 
probabilistic  model. 


there  are  fewer  targets.  In  fact,  when  the  number  of  targets 
grows  to  three  target  or  more  (Scenarios  C  and  D),  the 
localization  error  is  much  poorer  than  when  the  number  of 
targets  is  two  or  less.  The  CPF  exhibits  similar  problems  in 
Scenario  C  and  D  when  employing  the  disc  model  as  it  does 
for  the  other  two  scenarios.  Overall,  the  PHD  performs  very 
favorably  to  the  clairvoyant  filters  despite  the  fact  that  the 
PHD  is  actually  estimating  the  number  of  targets  in  the  scene 
during  each  snapshot  interval  when  the  disc  model  is  involved. 
Surprisingly,  the  clairvoyant  PHD  demonstrates  little  to  no 
advantage  over  the  PHD  for  any  of  the  sensing  model  cases 
over  all  scenarios. 

Figure  11  shows  the  average  estimated  number  of  targets 
versus  snapshot  time  for  the  PHD  filter  for  all  four  sens¬ 
ing  model  cases  over  all  four  scenarios.  For  the  matching 
probabilistic  case  (Case  1),  the  PHD  is  able  to  accurately 
estimate  the  number  of  targets.  When  the  PHD  still  uses  the 
probabilistic  model  but  the  sensor  behaviors  follow  the  disc 
model,  the  PHD  can  still  accurately  estimate  the  target  number 
as  long  as  the  targets  are  sufficiently  spaced.  When  the  targets 
are  nearby,  the  sensing  range  of  the  sensors  do  not  increase, 
and  the  size  of  the  cluster  of  sensors  detecting  the  sensors 
is  consistent  for  a  single  target  in  the  probabilistic  model. 
Therefore,  the  PHD  underestimates  the  number  of  targets. 
When  the  PHD  employs  the  disc  sensing  model,  the  ambiguity 
issues  related  to  the  disc  models  makes  it  impossible  for  the 
PHD  to  accurately  estimate  the  number  of  targets  over  all 
snapshots. 

B.  Effects  of  Innovative  Particles,  Motion  Model,  and  Initial¬ 
ization 

In  all  the  simulations  presented  up  to  this  point,  about 
9.09%  of  the  multiobject  particles  are  innovative  (J  =  100, 
Pex  =  1100).  Some  innovative  particles  are  necessary  to 
regulate  against  a  small  number  of  particles  and/or  a  simple 
motion  model.  In  the  next  set  of  simulations,  J  is  allowed  to 
vary  to  study  the  effect  of  the  number  of  innovative  particles  on 
performance.  Fig.  13  plots  the  performance  of  all  three  filters 
as  percentage  of  the  innovative  particle  when  Pex  =  1500. 
This  figure  represents  the  matching  probabilistic  sensing  case 
of  Scenario  A.  Similar  results  were  obtained  over  the  other 
scenarios  for  sensing  Case  1.  Note  that  when  J  =  1500, 
all  the  multiobject  particles  are  innovative,  and  the  motion 
model  is  ignored.  The  filters  are  simply  fusing  the  current 
set  of  binary  sensor  measurements  to  localize  the  targets. 
As  J  increases  from  zero  to  1500,  the  biters  are  simply 
putting  more  conbdence  on  the  current  measurements.  The 
bgure  demonstrates  that  about  10%  of  innovative  particles  are 
sufficient  to  accurately  estimate  the  number  of  targets  and 
their  locations.  The  performance  degrades  only  slightly  as  the 
innovative  percentage  grows  from  10%  to  90%.  The  local¬ 
ization  performance  falls  off  by  a  factor  of  two  when  using 
all  innovative  particles.  Thus,  the  tracking,  i.e.,  incorporation 
of  prior  measurements,  does  provide  for  some  performance 
advantages. 

The  good  performance  of  using  all  innovative  particles 
relative  to  using  none  implies  that  the  proposal  function  to 
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Fig.  11.  When  the  filters  implement  probabilistic  model,  estimated  number 
of  targets  using  PHD  for  (a)  Scenario  A  (b)  Scenario  B  (c)  Scenario  C  (d) 
Scenario  D.  The  legends  refer  to  the  actual  models  adopted  by  the  sensors. 
Note  that  for  the  matching  probabilistic  case,  the  average  RMS  error  of  the 
estimated  number  of  targets  over  the  whole  observation  period  is  0.34,  0.38, 
0.42  and  0.38  for  Scenario  A,  B,  C  and  D,  respectively. 


(d) 

Fig.  12.  When  the  filters  implement  disc  model,  estimated  number  of  targets 
using  PHD  for  (a)  Scenario  A  (b)  Scenario  B  (c)  Scenario  C  (d)  Scenario  D. 
The  legends  refer  to  the  actual  models  adopted  by  the  sensors. 
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(a) 


(b) 

Fig.  13.  For  Scenario  A  and  matching  probabilistic  case,  (a)  Average 
CA-OSPA  errors  (b)  estimated  number  of  targets,  versus  the  percentage  of 
innovative  particles  in  PHD 


obtain  innovative  particles  is  effective.  To  better  understand 
the  tradeoff  between  the  prediction  due  to  the  target  motion 
model  and  the  innovation  due  to  the  current  measurements,  in 
Fig.  14,  we  compare  PHDs  for  Scenario  1  and  Case  1  when 
one  uses  the  ground  truth,  the  random  samples  in  the  whole 
field,  and  the  random  samples  in  the  rectangular  area  bounding 
all  detecting  sensors  as  the  initial  PHD  particles,  respectively. 
Since  in  Scenario  A,  the  targets  are  moving  at  a  constant  veloc¬ 
ity,  the  ground  truth  initial  4-state  PHD  particles  are  the  tme 
target  positions  and  velocities.  For  the  random  initialization 
case,  the  initial  particles  are  random  samples  in  the  whole  field. 
For  the  detecting  sensor  initialization  case,  the  initial  particles 
are  random  samples  in  the  smallest  rectangular  bounding  box 
surrounding  the  detecting  sensors.  It  is  clear  that  the  ground 
truth  initialization  helps  the  most  when  there  are  no  innovative 
particles.  As  the  number  of  innovative  particles  increases,  the 
advantages  brought  by  the  ground  truth  initialization  become 
less  obvious.  When  using  the  detecting  initialization,  the  track 
filtering  employing  no  innovative  particles  is  as  effective  as 
employing  all  innovative  particle.  One  still  needs  to  use  around 
10%  innovative  particles  for  best  localization  performance. 
Note  that  the  smaller  range  for  the  y-axis  in  Fig.  14(a)  than 
in  Fig.  13(a)  shows  that  localization  performance  degrades 
slightly  as  the  number  of  innovative  particles  increases  from 
10%  to  90%.  Finally,  Fig.  14(b)  demonstrates  how  good 


(a) 


(b) 

Fig.  14.  For  Scenario  A  and  matching  probabilistic  case,  PHD  with  random 
initialization,  ground  truth  initialization,  detecting  sensor  initialization  (a). 
Average  CA-OSPA  errors  versus  the  percentage  of  innovative  particles;  (b). 
Average  CA-OSPA  errors  versus  time  when  the  innovative  particles  are  either 
zero  or  10  percent  during  tracking 


initialization  and/or  innovative  particles  enable  the  PHD  filter 
to  reach  a  steady  state  behavior. 

C.  Comparisons  with  ClusterTrack  in  1-D 

The  next  set  of  simulations  compare  the  PHD  against 
the  instantiation  of  ClusterTrack  as  described  in  [8].  Unlike 
the  PHD  method,  ClusterTrack  returns  particles  representing 
possible  tracks  where  a  particle  is  the  estimated  location  point 
of  time  at  each  snapshot  interval.  The  initial  iteration  of 
ClusterTrack  builds  up  the  particles  over  time  in  a  sequential 
fashion  so  that  at  time  t  —  T  a  particle  represents  the  estimated 
location  of  a  target  from  /  =  0  to  t  =  T  using  data  up  to 
time  T.  ClusterTrack  also  assumes  that  the  number  of  targets 
is  constant  and  estimates  the  target  number  as  the  maximum 
number  of  feasible  target  regions  over  the  entire  surveillance 
time  interval.  Then,  ClusterTrack  uses  multiple  iterations  to 
refine  the  particles  so  that  feasible  target  regions  are  covered 
by  the  collection  of  particles.  To  this  end,  ClusterTrack  does 
not  account  for  appearance  and  disappearance  of  targets.  Thus 
it  can  generate  more  particles  than  actual  targets,  and  its  ability 
to  capture  all  targets  depends  on  good  parameter  choices  which 
are  described  in  [8].  However,  it  is  admitted  in  [8]  that  there 
is  no  theoretical  guarantee  that  the  ClusterTrack  will  catch  all 
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the  targets.  Furthermore,  ClusterTrack  can  only  localize  all 
targets  after  multiple  iterations  which  is  operating  in  a  batch 
mode  where  measurements  from  future  times  are  influencing 
the  estimated  target  locations  at  earlier  points  in  time. 

In  our  simulations  for  the  ClusterTrack  method,  we  slightly 
modify  code  provided  by  the  authors  of  [8]  to  accommodate 
for  a  test  scenario  since  the  original  ClusterTrack  assumed 
the  number  of  targets  in  the  scene  remains  constant.  This 
first  test  reproduces  the  scenario  in  [8]:  sensors  are  uniformly 
distributed  in  1  -D,  the  separation  between  consecutive  sensors 
is  equal  to  the  sensing  range,  so  that  the  coverage  area  for 
the  adjacent  sensors  has  50  percent  overlap.  In  this  case,  the 
number  of  sensors  is  50,  and  the  disc  sensing  range  is  30,  and 
the  observation  range  is  1250.  Fig.  15  shows  that  estimated 
target  locations  from  the  PHD  filter  and  the  ClusterTrack  traces 
after  one  and  two  iterations.  The  PHD  and  ClusterTrack  (after 
two  iterations)  are  able  to  capture  all  targets.  However,  when 
ClusterTrack  is  behaving  as  a  filter  (and  not  a  smoother)  for 
the  first  iteration,  it  is  unable  to  capture  all  the  targets. 

In  the  second  comparison,  the  number  of  targets  changes 
during  the  surveillance  interval.  Figure  16  shows  that  the  PHD 
and  ClusterTrack  (after  four  iterations)  are  able  to  capture  all 
the  targets.  Again,  ClusterTrack  fails  to  localize  all  the  tracks 
after  the  first  iteration.  After  multiple  intervals,  the  traces 
in  ClusterTrack  will  jump  whenever  targets  are  born  or  die. 
ClusterTrack  forces  these  jumps  as  an  artifact  of  its  assumption 
that  the  number  of  targets  is  constant. 

Overall,  ClusterTrack  is  an  effective  multiple  target  tracker 
after  a  number  of  iterations,  and  the  localization  ability  of  the 
PHD  and  ClusterTrack  appear  to  be  comparable.  However, 
the  PHD  can  capture  all  targets  as  a  track  filter,  while 
ClusterTrack  requires  multiple  iterations  working  in  a  track 
smoothing  mode  to  capture  the  targets.  On  the  other  hand, 
ClusterTrack  does  provide  traces  of  the  tracks,  and  the  PHD 
only  provides  estimates  of  the  target  locations  at  each  snapshot 
interval.  Finally  the  method  to  build  and  refine  particles  in 
ClusterTrack  assumes  the  disc  model,  and  it  is  not  clear  how 
to  modify  ClusterTrack  to  consider  the  probabilistic  sensing 
model.  When  incorporating  the  probabilistic  sensing  models, 
the  PHD  provides  an  accurate  estimate  of  the  number  of  targets 
at  each  snapshot  interval  when  the  probabilistic  model  actually 
represents  the  behavior  of  the  sensors.  ClusterTrack  cannot 
provide  for  such  an  estimate. 

D.  Discussion 

The  simulations  demonstrate  that  effectiveness  of  the  pro¬ 
posed  PHD  method  for  multitarget  tracking  using  proximity 
sensors.  The  PHD  method  appears  to  be  more  robust  when 
employing  the  probabilistic  sensing  model  than  the  disc  model. 
Furthermore,  the  CPF  filter  can  break  down  when  employing 
the  disc  model.  It  appears  that  the  weight  update  using 
probabilistic  model  is  better  than  using  the  disc  model  because 
the  probabilistic  model  accounts  for  the  imperfect  multiobject 
particles  with  non-zero  weights.  Here  the  imperfect  multi¬ 
object  particles  refer  to  ones  where  some  subparticles  are 
good  but  the  others  are  far  from  any  ground  truth  target. 
Using  an  example  in  Fig.  17,  we  show  the  disadvantages 


(a) 


Iteration  =1 


(b) 

Iteration  =2 


(c) 

Fig.  15.  For  no  target  birth/death  cases,  (a),  estimated  target  positions 
using  PHD  vs.  time  (b)-(c).  at  each  iteration  of  ClusterTrack,  positions  of 
cluster  heads  vs.  time,  where  true  target  tracks  are  in  solid  lines,  and  dots  or 
dashed  lines  denote  the  estimates  using  either  PHD  or  ClusterTrack  method, 
respectively 
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Fig.  16.  For  target  birth/death  cases,  (a),  estimated  target  positions  using  PHD  vs.  time  (b)-(e).  at  each  iteration  of  ClusterTrack,  positions  of  cluster  heads 
vs.  time,  where  true  target  tracks  are  in  solid  lines,  and  dots  or  dashed  lines  denote  the  estimates  using  either  PHD  or  ClusterTrack  method,  respectively,  (f). 
the  number  of  targets  versus  time 


of  assuming  the  disc  model  in  the  measurement  update  for 
the  filters.  Following  the  measurement  update,  it  is  desirable 
for  the  weights  associated  to  imperfect  particles  to  be  greater 
than  zero  since  some  portion  of  its  subparticles  are  good. 
However  in  the  disc  model,  for  the  weights  to  be  nonzero,  all 
subparticles  associated  to  the  multiobject  have  to  lie  beyond 
the  sensing  radius  of  any  non-detecting  sensor,  and  there  must 
be  at  least  one  subparticle  within  the  detection  range  of  all 
detecting  sensors  (see  (5)  and  (13)).  For  example,  among 
the  five  sensors  in  Fig.  17,  three  sensors  centered  in  the 
shaded  circles  report  target  detections.  Suppose  there  are  three 
multiobject  particles,  each  denoted  by  the  stars,  pluses,  and 
diamonds,  respectively.  The  weights  for  the  star,  plus  and 
diamond  multitargets  are  0,  1,  and  0,  respectively,  when  using 
(5)  and  (13).  The  weight  for  the  star  particles  is  zero  because 
there  is  a  subparticle  that  lies  within  sensing  range  of  non¬ 


detecting  Sensor  4.  The  weight  for  the  diamond  particles  is 
also  zero  because  there  is  no  subparticle  that  lies  within  the 
sensing  range  of  sensor  3.  When  using  the  bounding  box  to 
generate  innovative  particles,  there  is  a  high  chance  that  none 
have  positive  weights  for  the  disc  model  when  the  targets  are 
far  apart. 

On  the  contrary,  the  weights  generated  by  the  probabilis¬ 
tic  model  for  the  imperfect  star  and  diamond  particles  are 
nonzero.  Overall,  the  probabilistic  model  does  a  better  job 
of  rating  the  overall  goodness  of  a  particle  as  the  weight 
is  not  binary.  Furthermore,  the  probabilistic  model  does  not 
prune  particles  that  contain  good  subparticles.  These  good 
subparticles  can  lead  to  good  solutions,  and  the  crossover 
feature  of  the  PHD  allows  good  subparticles  from  otherwise 
weak  particles  to  combine  into  strong  particles. 

When  the  number  of  targets  grows  to  three  or  more,  the 
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Fig.  17.  Illustration  of  weight  updates  using  the  disc  model  where  sen¬ 
sors!  1,2,3}  report  dt  =  1  and  sensors{4,5}  report  d[  =  0.  The  weights  for 
the  star,  plus  and  diamond  multitargets  are  0,  1,  and  0,  respectively,  when 
using  the  disc  model,  i.e.,  (5)  and  (13) 


Finally,  the  capability  for  the  PHD  methods  to  estimate 
the  number  of  targets  for  Case  1  is  not  boundless.  As  more 
targets  enter  in  the  scene  a  higher  percentage  of  sensors 
will  report  detections.  For  the  case  of  five  targets  (using  the 
probabilistic  sensing  model),  thirty  percent  of  the  sensors 
generate  detections.  When  the  number  of  targets  grows  to 
twenty,  all  hundred  sensors  will  normally  return  detections. 
In  this  extreme  case,  the  ability  to  localize  the  targets  com¬ 
pletely  disappears  as  any  configuration  of  twenty  targets  will 
be  sensed  by  all  sensors.  To  accurately  localize  targets,  a 
certain  number  of  non-reporting  sensors  are  required.  Perhaps, 
this  issue  also  affects  the  poorer  localization  performance  in 
Scenarios  C  and  D  when  targets  move  too  close  to  each  other. 
The  resolution  problem  does  not  appear  for  two  or  less  targets 
as  seen  in  these  simulations  and  in  the  results  of  [21]  that 
specifically  investigated  the  resolution  issue.  Future  work  can 
investigate  the  resolution  limits  of  the  PHD  filter  for  three  or 
more  targets. 


localization  accuracy  of  the  PHD  degrades  for  all  four  cases. 
However,  for  Case  1  (matching  probabilistic  model),  the  PHD 
is  able  to  accurately  track  the  number  of  targets.  It  is  interest¬ 
ing  to  note  that  for  Case  4,  the  filter  underestimates  the  number 
of  targets  for  Scenarios  C  and  D  once  the  number  of  targets 
grows  to  three  or  more  (see  Figure  11).  Since  the  PHD  filter 
is  employing  the  probabilistic  model,  the  underestimation  is 
caused  by  targets  that  are  too  close  to  each  other  to  create  large 
enough  clusters  of  detections  indicative  of  the  true  number  of 
targets.  As  a  result,  the  particles  of  the  PHD  will  coalesce 
into  a  smaller  number  of  clusters  than  targets  for  any  of  the 
four  cases,  and  this  will  create  problems  for  the  A:- means 
algorithm  to  determine  which  PHD  particles  should  belong  to 
which  targets.  Overall,  the  centroid  for  the  partitions  may  not 
be  reflective  of  the  ground  truth  targets.  Future  work  should 
investigate  whether  alternative  clustering  method  [29],  [30], 
[31]  can  improve  the  localization. 

In  all  simulations  in  Section  IV-A,  the  clairvoyant  PHD 
does  not  clearly  outperform  the  PHD.  This  was  unexpected 
given  that  the  clairvoyant  filter  always  forces  rv-k>  in  Step  3 
of  the  PHD  (see  Figure  2)  to  equal  the  ground  truth  number 
of  targets.  The  clairvoyant  PHD  does  slightly  outperform 
the  PHD  for  the  matching  probabilistic  sensing  case  for  all 
scenarios  except  Scenario  D,  when  four  or  more  targets  enter 
the  scene.  Inspection  of  Figure  11  indicates  that  the  PHD 
slightly  overestimates  the  number  of  targets  for  that  case, 
which  means  more  of  the  P*  PHD  particles  after  collapsing 
the  multiobject  density  (see  (23))  represent  good  solutions 
for  the  PHD  filter  than  for  the  clairvoyant  PHD  filter.  This 
larger  diversity  of  good  PHD  particles  might  translate  to 
better  localization.  In  fact,  forcing  n/1'1  to  equal  the  ground 
truth  target  number  does  constrain  the  search  for  good  mul¬ 
tiobject  particles  with  high  likelihood,  i.e.,  good  hypotheses 
that  explain  the  sensor  measurements,  and  may  offer  little  to 
no  advantages  when  the  sensor/target  configuration  is  such 
that  the  target  count  is  ambiguous.  Thus,  the  PHD  filter  can 
outperform  the  clairvoyant  PHD  filter  in  cases  where  the  disc 
model  is  considered. 


V.  Conclusions 

This  work  documents  our  investigation  of  multiple  target 
tracking  filters  in  proximity  sensor  networks.  To  this  end,  we 
developed  a  novel  formulation  of  the  PHD  filter  for  proximity 
sensors.  We  compared  the  PHD  against  two  clairvoyant  filters 
that  assume  the  number  of  targets  is  known  a  priori.  For 
measuring  the  localization  errors,  we  considered  a  measure 
that  attempts  to  separate  out  the  effects  of  cardinality  errors. 
The  simulations  show  the  advantages  of  the  crossover  feature 
in  the  PHD. 

Two  possible  measurement  models  were  considered:  the 
disc  and  probabilistic,  so  that  two  model  matching  cases 
and  two  model  mismatching  cases  were  considered.  While 
the  disc  model  has  wide  appeal,  the  probabilistic  model 
better  represents  the  phenomenology  that  dictates  how  sensors 
respond  to  target  signals.  The  simulations  demonstrate  that  the 
PHD  filter  would  be  unable  to  accurately  estimate  the  number 
of  targets  if  the  disc  model  was  assumed  in  the  measurement 
update.  This  is  the  case  because  the  hypothesis  that  explains 
the  pattern  of  detection  reports  is  not  unique.  On  the  other 
hand,  the  PHD  is  accurate  when  using  the  probabilistic  model 
in  the  measurement  update  because  of  the  uniqueness  of  the 
best  hypothesis  to  explain  the  pattern  of  detection  reports. 

The  unique  hypothesis  emerges  because  it  is  assumed  that 
the  radiating  power  of  each  target  is  a  fixed  and  known 
quantity.  In  reality,  one  would  expect  that  this  power  fluctuates 
from  target  to  target.  If  these  fluctuations  are  unbounded,  the 
probabilistic  model  becomes  completely  ambiguous,  and  it 
would  be  impossible  for  the  PHD  filter  to  accurately  estimate 
the  number  of  objects  in  the  scene.  However,  we  expect  that 
these  fluctuations  are  bounded.  Future  work  will  investigate 
how  allowing  for  bounded  spread  of  target  power  levels 
will  affect  the  performance  of  the  PHD  filter.  Other  future 
work  will  investigate  better  dynamical  models  to  improve  the 
tracking  performance  including  target  motion  and  birth/death 
models,  and  other  cluster  methods  for  closely  spaced  targets. 

Finally,  the  new  PHD  filter  in  this  paper  uses  a  profoundly 
different  measurement  model  than  the  traditional  PHD  filter. 
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When  the  sensor  measurements  are  statistically  independent 
conditioned  on  the  target  hypothesis  as  in  this  paper,  the 
multisensor  version  of  the  PHD  is  straightforward.  However, 
for  the  traditional  PHD  filter,  the  “exact”  multisensor  version  is 
very  cumbersome  to  implement  [32],  Usually,  approximations 
are  incorporated  in  the  multisensor  implementation  of  the 
traditional  PHD.  These  approximations  are  necessary  due 
to  the  constraints  that  a  measurement  corresponds  either  to 
one  unique  target  or  clutter.  Actually,  the  measurements  are 
usually  derived  from  raw  signals  through  signal  processing, 
e.g.,  direction  of  arrival  estimation.  As  an  alternative,  the 
new  PHD  method  could  easily  be  modified  to  consider  the 
likelihood  of  raw  signals  given  the  target  hypotheses.  The  new 
PHD  method  does  suffer  from  approximations  in  the  particle 
implementation  of  the  multiobject  target  density.  On  the  other 
hand,  approximations  of  the  PHD  for  sensor  fusion  also  lead  to 
performance  loss.  Future  work  can  investigate  whether  the  new 
PHD  method  operating  over  raw  signals  from  multiple  sensors 
has  utility  compared  to  multiple  sensor  implementation  of  the 
traditional  PHD  method. 
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